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We perform nonlinear general relativistic ideal magnetohydrodynamic simulations of poloidal mag- 
netic fields in rotating polytropic neutron stars. We have three primary goals: i) to understand the 
nature of magnetohydrodynamic instabilities inherent to poloidal magnetic fields in non-rotating and 
rotating neutron stars, ii) to explore the possible space of stable equilibrium configurations and iii) 
to understand gravitational wave emissions caused by the catastrophic reconfiguration of magnetic 
fields associated with giant magnetar flares. Our key physical contributions can be summarized as 
follows: i) gravitational waves from /-modes caused by magnetar flares are unlikely to be detected 
in the current or near-future generation of gravitational waves observatories, ii) gravitational waves 
from Alfven waves propagating inside the neutron star are more likely candidates, although this 
interpretation relies on the unknown damping time of these modes, iii) any magnetic held equilibria 
derived from our simulations are characterized as non-axisymmetric, with approximately 65% of 
their magnetic energy in the poloidal fleld, iv) rotation acts to separate the timescales of different 
instabilities in our system, with the varicose mode playing a more major role due to a delayed 
kink instability and v) despite the slowing growth rate of the kink mode, it is always present in 
our simulations, even for models where the rotational period is of the same order as the Alfven 
timescale. 

PACS numbers: 04.30.Db,04.40.Dg,95.30.Sf 



I. INTRODUCTION 



How can we probe magnetic fields in the cores of neu- 
tron stars? Despite their central role in multiple aspects 
of neutron star physics, this question has eluded sufficient 
resolution for nigh on five decades. Two independent sets 
of observations - neutron star spin down and thermal 
emissions - probe the magnetic field above the surface 
of the neutron star. However, these provide little infor- 
mation about the field lying inside the crust and core of 
the star. It is therefore left to first-principals modelling 
of neutron star interiors to garner these vital pieces of 
information. 

From a theoretical perspective, knowledge of the 
strength, topology and dynamics of magnetic fields in 
the core and crust of neutron stars is crucial to under- 
standing phenomena in garden-variety pulsars including, 
but not limited to; boundary conditions feeding magne- 
tospheric models that describe emission and neutron star 
spin-down [iH^, glitch dynamics |5| and igravitational 
wave emissions from internal fields [e.g. m]^ and mag- 
netic "mountains" [e.g. |^, [l^. Moreover, the class of 
neutron stars known as the magnetars, with surface field 
strengths oi B > 10^"* G, exhibit more exotic phenom- 
ena, whose existence and dynamics crucially depend on 
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the magnetic field. These include their unusually high 
surface temperatures [e.g [lll - [l3| (although see the recent 

article of Ho et al. Il3), the generation of magnetar fiares 
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1151 . Ila and their subsequent quasi-periodic oscilla- 
tions [e.g. [17h27| and free-precession [2g or lack thereof 
|29| . These reasons and more provide sufficient impe- 
tus to warrant first principals modelling of neutron stars 
emphasizing the effects of their magnetic fields. 

In two previous short articles we have begun to probe 
the interior dynamics of neutron stars strong magnetic 
fields utilising our newly developed three-dimensional, 
non-linear, general relativistic magnetohydrodynamics 
code HORIZON [ISl- The first of these [ly was a letter 
investigating a canonical series of non-rotating models 
with a single equation of state, in a first attempt to un- 
derstand poloidal magnetic field instabilities in general 
relativity. This letter explored the subsequent pseudo- 
equilibria states, broadening the possible solution space 
of stable equilibria available in neutron stars. A follow- 
ing short article jS^] utilised these models to calculate 
the gravitational wave emission from the complete re- 
construction of the magnetic field due to the magnetic 
field instability. This model was viewed as mimicking 
the behaviour of the interior of a neutron star imme- 
diately following a giant magnetar flare, and therefore 
gave estimates on the gravitational wave detectability of 
such a situation. In the present work we follow-up on 
those two short articles, providing significantly more de- 
tails on our numerical model, as well as confirming and 
extending our previous models using both different poly- 
tropic equations of state and also introducing rotation 



into the system. With this in mind, the present paper 
has three key, short-term motivations: i) to understand 
the nature of magnetohydrodynamic (MHD) instabihties 
inherent to poloidal magnetic fields in both non-rotating 
and rotating relativistic neutron stars, ii) to substantiahy 
extend the full set of MHD equilibria derived as steady- 
state solutions and iii) to robustly extend the results re- 
garding gravitational wave emission following giant flares 
to include different polytropic equations of state, and 
hence include a relationship between the gravitational 
wave strain, magnetic field, radius and mass of the star. 

The study of axisymmetric MHD equilibria relevant 
for neutron stars has a rich heritage, established with the 
early work of Chandrasekhar and Fermi [33] , Monaghan 
[3J|, Roxburgh |35| and Parker |36|. It was soon realised 
that these idealised analytic fields are prone to various 
kinds of instabilities including, but not limited to, Tayler 
"kink" instabilities |37h41| and the Flowers and Ruder- 
man (421 instability. In their various guises, these insta- 
bilities indicate that purely poloidal and purely toroidal 
fields are dynamically unstable. Poloidal fields are unsta- 
ble in the regions where its field lines are closed within 
the star, and it is therefore widely believed that threading 
a toroidal component through this region of the poloidal 
field acts to stabilise the field [39, 43]. Such twisted-torus 
configurations have been studied in considerable detail in 
semi-analytic calculations [y, |44| - |47| and also derived as 
equilibrium configurations from global MHD simulations 
p8, 49]. However, the stability of such configurations in 
barotropic stars has recently been questioned [50]. 

In our recent paper beginning with axisymmetric, 
purely poloidal initial conditions [3l|, our evolutions 
characteristically developed non-axisymmetries during 
the early phases of the non-linear development of the 
kink instability. All of our subsequent equilibria where 
therefore non-axisymmetric, although they retained cer- 
tain properties of the twisted-torus configurations. The 
development of these non-axisymmetries were consistent 
with the evolutions of Braithwaite 51], who found that 
the evolution of fields to non-axisymmetric or axisym- 
metric states depended on the radius of the initial neutral 
line. 

Rotation has a marked effect on the development of 
various MHD instabilities. From an early stage, Frieman 
and Rotenberg [52] understood that rigid-body rotation 
has an effect on the stability of hydromagnetic equilibria 
only when the velocity of the fiuid fiow is of the same or- 
der, or greater than the Alfven velocity. The first works 
on the global stability of purely poloidal, rotating fields 
was that of Geppert and Rheinhardt [53] , who found that 
stars rotating with sufficiently high rotational velocities 
have suppressed instabilities. Lander and Jones |54'|, on 
the other hand, utilised their linear code to conclude that 
some, but not all modes were stabilised by the presence 
of rotation. In the present article, we perform the first 
general relativistic simulations of rotating poloidal fields, 
concluding that instability timescales are slowed by the 
rotation, but that the instabilities are not completely 



suppressed. We do note, however, that we are only be- 
ginning to approach the regime of fast rotation, whereby 
the fluid-flow velocity is of the same order as the Alfven 
velocity. We leave open the distinct possibility that faster 
rotation would suppress these instabilities. 

Throughout the article we allow Greek indices to range 
... 3 and Latin indices 1 ... 3. The paper is set out as 
follows: in section |lT] we outline our numerical method, 
focussing on the equations of general relativistic magne- 
tohydrodynamics in III Al and details of our specific code 
implementation in IIIBl In section IIIII we look at a non- 
rotating fiducial model, concentrating on the magnetic 
field instability in IIIIBl gravitational wave emissions in 
IllICi quasi-equilibrium end-states in section IIIIDI and 
the effect of different magnetic field strengths in IIII El In 
section HVl we study the relationship between the radius 
and mass of the star, magnetic field strength and grav- 
itational wave emissions, deriving for the first time an 
empirical relation between these four variables. In this 
section we also look at the detectability of such gravi- 
tational waves in current and future gravitational wave 
detectors such as Advanced LIGO and the proposed Ein- 
stein Telescope respectively, concentrating on both the 
/-mode emission and also lower frequency Alfven modes. 
In section |V] we look at the effect rotation has on the 
varicose and kink instabilities. We conclude in section 
lYll 



II. NUMERICAL MODEL 

We are studying magnetised neutron stars through the 
time evolution of the ideal MHD equations in general rel- 
ativity utilising the HORIZON code [30l - [3^ . In Zink [s^] 
and Lasky et al. [3l| we have presented brief outlines of 
our numerical method, however it is worth elaborating 
on this in more detail. Throughout this section we as- 
sume geometrised units such that G = c = 1, although 
we retain specific units for the remainder of the article 
following this section. 



A. General Relativistic Magnetohydrodynamics 

To express the equations of General Relativistic Mag- 
netohydrodynamics (GRMHD) in a form appropriate for 
numerical integration, we follow closely the formalism 
outlined in Gammie et al. |55| . To this end, we define 
a four-vector that is orthogonal to hypersurfaces of con- 
stant time, t, which has components "-^ = ^ (l, —/?')■ 
Here, a is the lapse function and /3* are the spatial 
components of the shift vector. The line element of 
such a spacetime can then be expressed in coordinates 
x'^ = (t, x*) as 

ds^ = - (a2 _ /3„/3") dt^ + 2Padtdx'' -f -/atdx^dx'', (1) 

where 7^1/ = 5^^ + n^n,y is the induced three-metric. 



We denote the four- velocity comoving with the fluid as 
u^, allowing a definition of the three- velocity of the fluid 
as measured by the observer moving along n^; 






(2) 



where W 



-UaTL 



is the relative Lorentz factor between 



the two observers and hij is the induced three-metric on 
the constant t hypersurfaces. 

The electromagnetic field is now defined by the anti- 
symmetric Faraday tensor, F^j,^, with the four Maxwell 



equations given by 



Vr„,K 



[/i-f^i/o-] 



=0, 



V"F^„ ^inj^. 



(3) 

(4) 



Here, JT'^ is the electric four-current and square brack- 
ets denote anti-symmetrization of the indices. The four- 
current can be decomposed into components in and or- 
thogonal to the fluid four-velocity which, together with 
Ohm's law, can be expressed as 



Jt" = eu" + aPf^^Ua, 



(5) 



where e is the proper charge density and a is the electric 
conductivity as measured in the comoving frame. 

The electric and magnetic fields for any observer can 
be evaluated by appropriately contracting the Faraday 
tensor with the observers four velocity. For example, 
the electric and magnetic field vectors as observed in the 
frame orthogonal to hypersurfaces of constant t can be 
expressed respectively as 






AiaP'i 



-^Q^^7l 



(6) 
(7) 



where *F^^^ is the Hodge dual of the Faraday tensor and 
g/ii/o-T jg |-j^g Levi-Civita alternating pseudo-tensor. 

Moreover, we define the magnetic field observed in the 
comoving frame of the fiuid as 



6^=*F^"u„. 



(8) 



Throughout this article we assume our fiuid to be per- 
fectly conducting (i.e. we are working in the ideal MHD 
approximation) , implying we take the limit of cr — > oo . In 
order to keep the current finite, this implies -F^'^Mq — 0, 
implying the electric field observed by the comoving ob- 
server vanishes. This property allows the electric field in 
any frame to be expressed in terms of the magnetic field 
and the relevant four- velocities, implying the electric field 
no longer enters the calculations. 

The stress-energy tensor is expressed in terms of the 
fiuid (assumed herein to be a perfect fluid) plus electro- 
magnetic components, which are respectively given by 



Tuv"^ = phu^Uy+pg^^, 



rpEM 
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(9) 
(10) 



Here, p is the energy- density of the fluid, h = 1 + e+p/ p 
the specific enthalpy, with e being the specific internal 
energy and p the isotropic pressure. We assume a poly- 
tropic equation of state (EoS) for our neutron star mod- 
els, such that p = Kp^ , discussing this in more detail 
below. 

With the above stress-energy tensors, the conservation 
law VqT"^ = 0, can be expressed in a coordinate basis 
as 



d 



TO iV^Tf,') = y^T^^r"^^, 



(11) 



where F^j^o- are the Christoffel symbols. Moreover, the 
spatial and temporal components of the induction equa- 
tion can also be expressed respectively as 



[V^{Vu'~b'u^)]=0, (12) 
d 



r—gB^) =0. (13) 



B. The HORIZON Code 

Horizon is a GPU based numerical code borne out of 
the CPU, general relativistic hydrodynamics code thor 
|56l |57| . As such, horizon solves the equations of 
GRMHD outlined in section Hi Al Throughout this arti- 
cle we generate initial conditions using the lorene spec- 
tral code^ , which produces self-consistent solutions of the 
coupled Einstein-Maxwell field equations in ideal MHD 
|58| . The LORENE solver finds only solutions with purely 
poloidal magnetic fields including rotation, and we there- 
fore restrict our attention to these models in the present 
article^ . 

Following the generation of initial conditions, we map 
the spectral grid to a regular, three-dimensional Carte- 
sian mesh. Throughout the present article we do not add 
any artificial perturbation to the system, relying instead 
on noise created through the mapping process to add a 
pseudo-random perturbation to the system. Our Carte- 
sian grid typically contains 120^ grid points. We have 
extensively tested this grid resolution, with some details 
given below. We impose a low-density, artificial atmo- 
sphere in the region of our spacetime exterior to the star. 



http://www.lorene.obspirL.fr/ 

It is worth noting that alternative methods exist for generating 
initial conditions, for example the numerical codes of Bucciantini 
and Del Zanna _59i known as XNS] or Kiuchi et al. 60], which 
both find axisymmetric equilibria for stars with purely toroidal 
fields. Alternatively, one can use the RNS code 61] for solving 
the equations of general relativistic hydrodynamic equilibrium, 
subsequently imposing a magnetic field as a perturbation. This 
latter approach, while initially providing a small constraint vio- 
lation, has the advantage that arbitrary magnetic field topologies 
can be imposed on the system. 



That is, for all grid cells for which the density drops be- 
low some critical value, we impose p = Patmj which for 
the evolutions presented herein is Patm = 10^^ in units 
where c = G = Mq = 1. We subsequently allow for 
the full evolution of the magnetic field in this region (see 
below for details of the evolution). This is in contrast 
to many GRMHD numerical studies where, for exam- 
ple, the magnetic field is confined to the interior [e.g., 
|62| . |63| . the magnetic field is prohibited from evolving 
in the exterior [6Jj or where ad hoc magnetic diffusivity 
terms are added to the induction equation [6^, [6a] . Our 
method still does not treat the exterior of the star cor- 
rectly (one would require significantly more complicated 
magnetospheric physics such as radiation transport and 
force- free magnetic fields), however it does allow for a 
free evolution of the magnetic field at the surface of the 
star which is important for the dynamical evolution. 

The outer boundary of our star is located approxi- 
mately 1.4 times the stellar radius (at the closest point), 
and there we adopt Dirichlet boundary conditions for the 
evolution of the magnetic field. Dirichlet boundary con- 
ditions are restrictive in the sense that they do not al- 
low the magnetic field to evolve at the outer boundary. 
However, we persist with these conditions as they are 
found to be the most stable for our numerical simula- 
tions. In particular, we have attempted both Neumann 
boundary conditions as well as linear extrapolation tech- 
niques, both of which eventually introduce numerical in- 
stabilities into the atmosphere of the star. For a limited 
subset of these simulations we were able to evolve for long 
enough to track the magnetic field instability (see section 
imp , exhibiting minimal difference in growth timescale 
and topological behaviour of the instability between dif- 
ferent boundary conditions. We are unable to rigorously 
test our equilibrium configurations against these bound- 
ary conditions due to the numerical instabilites, however 
we do note that a majority of the magnetic field evolution 
is driven by the interior dynamics of the star, with the ex- 
terior magnetic field evolving very little. We discuss this 
in more detail in section IIII Dl It is also worth noting 
that we have rigorously tested our evolutions against the 
location of the outer boundary and found no discernible 
difference with either the nature of the instability or the 
end-state of the magnetic field configuration. 

Throughout the present article we adopt the Cowling 
approximation, such that the spacetime metric is held 
fixed throughout the evolution. This has the obvious ben- 
efit of significantly reducing computational costs. More- 
over, throughout the article we shall be looking at recon- 
figurations of the magnetic field, which generally occurs 
on small-scales around the stellar interior. In Lasky et al. 
[31j we discussed how these instabilities act on equipoten- 
tial surfaces, implying Cowling is a good approximation. 
We do note however, that the adoption of the Cowling 
approximation does not allow us to study the damping of 
various oscillation modes due to gravitational wave emis- 
sion. Particularly for the /-mode, this is anticipated to 
be the main source of damping, acting on timescales of 



order 0.1s [67[. Other, lower frequency modes however 
may last significantly longer [6^, a point discussed in 
more detail below. 

The initial conditions are given in terms of the prim- 
itive variables, and must therefore be converted into 
the conservative variables utilised in the evolution code. 
These relations are algebraic, and can therefore trivially 
be evaluated in each cell. The temporal evolution then 
requires the computation of the cell face flux vectors, 
which are obtained by solving a local Riemann problem 
across each cell (see Ref.j30| for details). This method 
solves for the conserved variables at each time step, how- 
ever the source terms require knowledge of the primitive 
variables at each step. Transforming the conserved vari- 
ables to the primitive variables is not, in general, an al- 
gebraic process. To this end we adopt a one-dimensional 
polytropic recovery scheme [69|, employing a Newton- 
Raphson method to solve the non-linear algebraic equa- 
tions. As a final note, we adopt hyperbolic divergence 
cleaning according to the prescription outlined in Ander- 
son et al. [7(11 , which maintains a divergence free magnetic 
field throughout the evolution. 



C. Shock Tube Code Tests 

Prior to detailing results, we shall spend some time 
exploring two standard shock tube code tests which 
give confidence in the numerical accuracy of our results. 
These are the Sod [7l| and Balsara [72| tests, which were 
first presented for the horizon code in Zink |3g |. 

The Sod [7l| test only includes the hydrodynamic por- 
tion of the code (i.e. in the absence of magnetic fields) . In 
particular, a Riemann problem is prepared for an ideal 
gas with r = 5/3 and initial states pi = 1, Pl = 1, 
< = and PR = 0.125, Pr = 0.1, u^ = 0. This prob- 
lem is prepared on a full 120'^ grid, and simulated up to 
t — 0.8. A reference solution has further been created 
using the well-tested thor code 56^, _57J with a very fine 
grid. 

In figure [T] we plot the density (top panel) and the x 
component of the velocity profiles at the end of the evo- 
lution for both single and double precision simulations, 
as well as for the reference solution. The dissipative na- 
ture of the numerical method can be seen with the subtle 
difference between the horizon runs and the finer-grid 
THOR simulation. The largest difference between the sin- 
gle and double precision results are \6p\ sa 5 x 10~^ and 
l^u^l w 10~*, which appear exactly at the location of the 
shock. In the smooth parts of the fiuid flow, these differ- 
ences are typically \Sp\ < 10^^ and \6u^\ < 10~^. These 
errors are at or lower than the level of the discretisation 
error of the numerical scheme. 

The Balsara [73 test includes the magnetic fleld, which 
contains the same initial conditions as for the Sod test, 
with the additional conditions of B]^ = (0.51, 0) and 
^R — (0-^! ~1' 0)- I^ figure [2] we plot the y component 
of the magnetic field vector at the end of the evolution 
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FIG. 1: Density profile (top panel) and x component of the 
velocity field («^; bottom panel) for the Sod test. The red 
line with open circles and the blue thin line are the hori- 
zon evolutions using single and double precision respectively. 
The reference solution (thick black line) was obtained with a 
very high resolution simulation by the well-tested thor code. 
These plots show the end of the simulation after t — 0.8. 
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FIG. 2: Profile of the y component of the magnetic field, B^ , 
for the Balsara test after t — 0.8 showing both single and 
double precision accuracies. The reference solution is again 
obtained using a very high resolution grid in the thor code. 



for single and double precision and also for the reference 
solution. The magnetic field evolution is not seen to be 
affected by the use of single precision accuracy, with the 
absolute difference being \SBy\ < 10^^. Similar values 
hold for other evolved variables. 

Details of tests on rapidly rotating, non-magnetised 
neutron stars, as well as tests involving toroidal magnetic 
fields in non-rotating magnetised stars, can be found in 



Zink [3(y\. These include recovering the correct spectral 
frequencies (as compared with accurate, linear codes) of 
various fluid modes for rapidly rotating stars. We omit 
the details of these tests in the present article as we shall 
also be recovering the frequencies of various oscillation 
modes in our fiducial models and, in particular, compar- 
ing these with the results from accurate linear codes (see 
section UlI A[) . Moreover, we detail the remainder of our 
numerical tests performed for our specific fiducial mod- 
els in the following sections alongside the corresponding 
results. 



III. FIDUCIAL SIMULATION 

In this section we review the key results of Lasky et al. 
[3l| and Zink et al. [32], providing significantly more de- 
tails than published in those two works. In particular, in 
order to discuss the magnetic field kink and varicose in- 
stabilities, as well as the equilibrium configurations pre- 
sented in Lasky et al. [3l|, we utilise a fiducial model 
with poly tropic EoS with ii' = 100, F = 2. Such a 
model has gravitational mass of 1.3 M© and equatorial 
radius R — 12.6 km. Moreover, this fiducial model has 
an average internal magnetic field strength (i.e. where 
P > Patm) of Bi5 = 13, where B15 = B/10^^ G. To 
make contact with neutron star observations we express 
the magnetic field in terms of the surface field strength 
evaluated at the pole of the star, which for this fiducial 
model is B15 = 8.8. Finally, the characteristic Alfven 
timescale, defined according to 



TA = 



2R^A^ 



(B) 



(14) 



where (. . .) represents a volume weighted average, evalu- 
ates to TA = 5.0 ms for our fiducial models. 

In figure[3]we plot the evolution of the rest mass for our 
fiducial simulation. One can see an initial loss of mass 
in the first few milliseconds of the simulation. This cor- 
responds to errors associated with mapping the lorene 
spectral grid to our Cartiesian grid. This initial mass loss 
soon reaches a new equilibrium, at which point the mass 
evolves almost unchanged until approximately 75 ms into 
the simulation. At this point, which corresponds to the 
nonlinear saturation of the magnetic field instability (see 
section IlIIBp . the mass begins to slowly increase. Af- 
ter more than t = 360 ms we see a total change in the 
rest mass of less than 0.15%. We note that, even for 
our models with the strongest magnetic fields of more 
than 10^^ G in the center of the star, we see less than a 
0.3 % change in the total rest mass over many hundreds 
of milliseconds. 

In figure [H we present evolutions of the rest mass for 
different resolution simulations, using the same fiducial 
model presented in figure [31 In particular we show four 
simulations with 90'^, 120'^ (our canonical model from fig- 
ure [3]), 150"^ and 180"^ grid-points. One can see here that 
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FIG. 3; Rest mass as a function of time for a non-rotating 
polytropic model with K = 100, F = 2 and central magnetic 
field of 1.0 X 10^^ G. An initial rearrangement of the rest mass 
is seen in the first milliseconds. This quantity is conserved on 
the level of ~ 0.3%. 
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FIG. 4: Evolution of rest mass for model presented in figure 
[3] with 120"^ grid-points at single precision (SP; black thick 
line), 120 grid-points at double precision (DP; red dotted 
line) and 180'' grid-points at SP (dashed blue line). 



the lower resolution simulations have a larger mapping 
error at the beginning of the simulation, however over 
the remaining 70 ms shown in this plot there is little de- 
viation in the rest mass, independent of resolution. 



A. Mode Analysis 

The benefits of understanding pressure related modes 
in our fiducial model are two-fold. Firstly, they provide 
an excellent code check that we are reproducing known 
results from the literature, and secondly they will aid in 
our understanding of the gravitational wave emission (see 
sections ImC] and |IV|. 

In figure [5] we plot the evolution of the central energy- 
density, pc as a function of time. As mentioned above, 
our simulations are only perturbed through the mapping 
process between the lorene spectral grid and our Carte- 
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FIG. 5: Evolution of the central energy-density, pc, for our 
canonical model. The full evolution lasts over 360 ms, while 
the inset in the plot shows only the first 15 ms of the evolution. 
Note that the central energy- density is given in natural units 
such that c — G — Mq = 1. The large perturbation seen 
in the initial evolution is due to the mapping between the 
LORENE spectral code and our Cartesian grid. 



sian grid on which the evolutions are performed. In the 
initial moments of the simulation, i.e. for approximately 
the first 15 ms, this induces a perturbation in the central 
energy-density that can be seen clearly in the inset of fig- 
ure [51 This initial perturbation is damped by numerical 
viscosity such that is it no longer visible approximately 
15 to 20 ms into the simulation. The central density then 
evolves with little variation until approximately 50 ms, 
at which point the magnetic field instability is seen to 
disrupt the density in the middle of the star. The subse- 
quent evolution of the central density also contains sig- 
nificant motion due to the additional kinetic energy in 
the star provided by the magnetic field instability. We 
discuss this in significantly more detail in the following 
section. 

In figure [S] we show the Fourier transform of the cen- 
tral density, pc, for the first 15 ms of the evolution (black 
line) and for the entire 360 ms (blue line). The two curves 
are scaled such that the sizes of the first peak are equiv- 
alent. The black curve provides a significantly cleaner 
signal than the blue curve due to the initial signal shown 
in the inset of figure \5\ The four strong peaks seen at 
approximately 2530, 4280, 5990 and 7750 Hz represent 
the purely radial, £ = 0, F-mode and its three lowest 
overtones respectively. 

We compare this _F-mode frequency with that found 
in the literature. In particular, Gaertig and Kokkotas 
[TSiI utilised a linear numerical code to study stellar os- 
cillations oi r — 2, K — 100 polytropes, however with 
different central densities to those used herein. They used 
stellar models with gravitational mass M = 1.4 M©, and 
equatorial radius R — 14.15 km, for which the F-mode 
was found to have a frequency of 2679 Hz. Given that the 
mode frequency scales with the compactness of the star, 
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FIG. 6: Fourier transform of the central energy-density, pc, for 
our canonical model. The black line is a Fourier transform of 
the first 15 ms of the evolution, which has the greatest pertur- 
bation, and hence cleanest signal for the pressure modes. The 
blue line is for the entire evolution lasting more than 350 ms. 
Here we see the fundamental radial f -mode at ~ 2530 Hz and 
its overtones. 



M / R, scaling their frequency to our stellar model gives a 
corresponding i<"-mode frequency of 2547 Hz, which can 
be compared to our derived result of 2530 Hz. This small 
difference in frequency is within our error margin given 
the calculation however we further note that the Gaertig 
and Kokkotas [73| calculations did not include magnetic 
fields as was done herein. 

As expected, the central energy-density evolution 
shows no hint of any non-radial pressure modes. To look 
at these modes we track various quantities at different 
radii throughout the star. In particular, for this, and 
also for the purpose of tracking the magnetic field insta- 
bility (^section FHIBp , we evaluate a Fourier decomposition 
of various physical quantities on a ring in the equatorial 
plane. That is, we compute complex weighted averages 
[e.g.lzl 



Cnril) 



1 

2^ 



/ (nj, 0, z = 0) e""'^c 



(15) 



where w — \/ x^ + y^ = const, lies in the initial equa- 
torial plane of the magnetic field. To track fluid modes 
we look at the quantity Cm (p)- In figure [7] we plot the 
Fourier transform of the real part of Co (p), extracted 
at ci7 = 0.75tx7j,, where Wi, is the equatorial stellar ra- 
dius. Again, the black line represents the Fourier trans- 
form of only the first 15 ms of the evolution, while the 
blue line is the Fourier transform of the entire evolu- 
tion. In this signal we see considerable noise in the lower 
frequency band of the spectrum, however we also see a 
considerable number of distinct modes. In particular, the 
first mode located at approximately 1750 Hz is the fun- 
damental, non-radial 1 = 2 /-mode. The second peak is 
again the F-mode seen in figure [51 and the final mode at 



FIG. 7: Fourier transform of the real part of Co (/o) as defined 
in equation (|15|l and extracted at three-quarters of the stellar 
radius. The black line is a Fourier transform of the first 15 ms 
of the evolution, which has the greatest perturbation, and 
hence cleanest signal for the pressure modes. The blue line is 
for the entire evolution lasting more than 360 ms. The fun- 
damental 1 = 2 non-radial /-mode is seen at ~ 1750 Hz, the 
F-mode at ~ 2530 Hz and the 1 = 2, pi-mode at ~ 4250 Hz. 



^ 4250 Hz is the pi-mode. 

Taking literature values from Gaertig and Kokkotas 
[73.] for the £ = 2, /-mode, we find a re-scaled value of the 
frequency for their models to be approximately 1797 Hz, 
which is again in rough agreement with the 1750 Hz fre- 
quency found herein. 



B. Instability 

A broad appreciation of the magnetic field instability 
can be gained from looking at the time evolution of the 
change in magnetic energy, plotted in figure [S] Here, 
AE/Eq = (E — Eq) /Eq, where E is the total magnetic 
energy and Eq is the magnetic energy a.t t = 0. We see 
an initial rearrangement of the magnetic field that settles 
after approximately t = 10ms = 2ta- This is associated 
with mapping errors present in the rest mass and energy 
density evolutions, however this takes longer to settle to 
an equilibrium as the timescale for the magnetic energy is 
the Alfven crossing time rather than the sound crossing 
time. 

Following the initial rearrangement, the magnetic en- 
ergy remains constant until approximately 50ms, at 
which point a sharp decrease is seen in the energy which 
is attributed to the kink instability (see below). This loss 
of magnetic energy is therefore associated with a conver- 
sion to kinetic energy in the system. 

This magnetic field instability can be studied in more 
detail by again computing complex weighted averages of 
the Fourier decomposition presented in equation (J15l) . 
This time we choose the function / to be the (f> com- 
ponent of the magnetic field, B^, although we note that 
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FIG. 8: Evolution of the change in total magnetic energy, 
AE/Eo, for our fiducial model. 
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FIG. 9: Evolution of Cm (B^) for our fiducial model with 
average magnetic field strength of Bis = 13, which has an 
Alfven crossing time of 5.0 ms and surface polar magnetic 
field of _Bi5 — 8.8. The arrows represent the times of the 
three-dimensional snapshots shown in figure [TO] and the two- 
dimensional equatorial slices in figure [12] 



other components of the magnetic field, fluid velocity and 
energy-density also exhibit the instability, however the 
signal is generally not as clean. 

In figure [5] we plot the evolution of the m — 1 , . . . , 4 
modes for Cm (B^) for our fiducial model. These quanti- 
ties are measured at tu = 0.6-07^, where tu* is the stellar 
radius. It is worth noting that the neutral line of these 
simulations (i.e. the ring around the equatorial plane 
where B = 0) is located at approximately two-thirds of 

In figures [TU] we present three-dimensional plots of 
the magnetic field at various instances throughout the 
evolution. The blue volume rendering in these figures 
is an isopycnic surface of p = 0.37pc, where pc is the 
central rest-mass density. This value was chosen as it 
lies approximately 50% of the radius of the star, and 
therefore provides contrast to the magnetic field lines. 
Two sets of magnetic field lines are plotted, both of 
which are seeded on the equatorial plane and traced both 



in the positive and negative B direction. The red field 
lines are seeded close to the neutral line in the equatorial 
plane, whereas the black field lines have been seeded 
interior to the neutral line. It is worth noting that the 
domain of our Cartesian grid is larger than that plotted 
in figure [TU] and that the field lines are truncated at the 
surface of the star for clarity in the images. A movie of 
the fiducial simulation lasting 400 ms can be viewed at 



http : //www . tat . physik . uni-tuebingen . d.e/~tat/grmhd/ 



Figure [TUk shows the initial data as imported from the 
LORENE spectral code, as described in section [IT] Figure 
[TUb shows the evolution of the system after t — 25 ms = 
5 TA ■ One can clearly see the onset of the 'sausage' or 
'varicose' mode [40] in this snapshot. The sausage mode 
involves a change in the cross-sectional area of a fiux tube 
around the neutral line. It is clearly visible that this 
has developed most strongly in the m = 4 mode, which 
is confirmed through a comparison with figure [9] which 
shows this mode clearly dominating over the others in 
this early phase of the evolution. 

The m = 4 varicose mode in these non-rotating simu- 
lations is the result of a transient excitation attributed to 
the Cartesian grid. We have verified that this transient 
reduces with increased grid resolution, however note that 
the presence of the varicose mode is an inherent charac- 
teristic of the system. The varicose mode is discussed in 
more detail in Lasky et al. 3l[ for non-rotating simula- 
tions, and is also discussed in greater detail in section \V\ 
of the present paper in the case of rotating models. 

The varicose mode visually dominates our simulations 
for almost ten Alfven crossing times before the 'kink' 
instability appears and begins to dominate the system. 
This is presented in figure [TUb . which is taken at t = 
50ms = 10 TA. One can see that the kink instability is 
acting orthogonal to the gravitational field [40] , with the 
presence of the varicose mode still clearly visible. This 
mode will have been excited from the beginning of the 
simulation (as seen clearly in figure [U]), however it is only 
at this point that the exponential growth has reached 
a point at which it is visually obvious. This therefore 
represents the non-linear development of the instability 
where the change in the field structure is of similar order 
to the background field. 

As discussed in Lasky et al. [3l|, the modal analysis 
presented in figure [9] does not distinguish between the 
varicose and kink modes, implying exponential growth 
does not indicate an instability in one or the other mode. 
In figure [11] we plot the ratio of magnetic energy in the 
poloidal field, Ep, to the total magnetic energy, E, as a 
function of time for our fiducial simulation. The initial 
configuration imported from LORENE is purely poloidal, 
i.e. Ep/E — 1. As discussed, the varicose mode is vi- 
sually apparent in the simulations for the first ~ 40 ms, 
although we can see from figure [TT] that this mode has no 
effect on the poloidal energy ratio. The non-linear devel- 
opment of the kink instability however, causes the Ep/E 




FIG. 10: Time evolution ol fiducial model with average magnetic field Bis = 13, corresponding to an Alfven wave crossing time 
of 5ms. The figures are a.) t — ms, h) t — 25 ms, c) t — 50 ms and d) t = 195 ms. To more clearly visualise the instability, the 
red field lines are seeded on the equatorial plane close to the neutral line, and the black field lines are seeded on the equatorial 
plane interior to the neutral line. The volume rendering is an isopycnic surface at 37% of the central rest-mass density, shown 
to provide contrast with the field lines. 



ratio to significantly change. In other words, it is only 
the kink mode that introduces a non-zero toroidal com- 
ponent of the field. This combination of modal analysis 
and magnetic energy ratios will become important for 
understanding the evolution of the rotating simulations 
in section IVl 

Finally, in figure llOd we show a typical quasi- 
equilibrium state of the simulation, which, in this case, is 
t = 195 ms = 39 ta into the simulation. We discuss this 



in more detail in section fill Dl 

Figures [To] can be somewhat misleading in an interpre- 
tation of the neutral line of the star. For this reason, 
in figures [T^ wc plot the absolute value of the magnetic 
field on an equatorial slice for the same four timesteps 
plotted in figure [TUl In particular, figure fT2k shows the 
initial conditions where one can see the neutral line at 
a constant radius of approximately 8.5 km (compared to 
the surface of the star at 12.6 km). Figure [T^ represents 
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FIG. 11: Ratio of poloidal-to-total magnetic field energy, 
Ep/E, as a function of time for the simulation presented in 
figures l9l and [TOl The initial configuration is purely poloidal 
{Ep/E = 1), with the kink instability shown in figure [TOb 
causing a large conversion of energy between poloidal and 
toroidal components. The instability saturates at approxi- 
mately Ep/E ~ 0.65. 



t = 25 ms = 5ta where the to = 4 varicose mode has 
come to visually dominate the system. One can see here 
that the size of the flux tube around the neutral line has 
remained constant, however the m — A excitation dis- 
cussed above has caused the shape of the neutral line to 
change. Figure [T^t is at i = 50 ms = 10 ta, at which 
point the kink and varicose modes can both be clearly 
seen. In the kink mode, the flux tube around the neutral 
line has kinked above and below the equatorial plane, 
which can be seen from the oscillatory nature of the field 
minima. Again, figure (T^Ji is after t — 195 ms = 39 ta 
and is discussed in more detail in section ITlID I 



Gravitational Wave Emission 



In Zink et al. (32| we utilised these simulations to place 
estimates on the gravitational wave emission due to an 
internal rearrangement of a magnetic field. That is, we 
used these simulations to mimic the aftermath of a giant 
magnetar flare that acts to rearrange the internal mag- 
netic field configuration (see also Ref. [66|). We calculate 
the gravitational wave strain from our simulations utilis- 
ing the quadrupole formula [e.g..l75| 



strain, at each time-step we evaluate lij, which can be 
expressed in terms of spatial quantities utilising the con- 
tinuity equation. This significantly decreases differenti- 
ation errors due to relatively large time-steps as we are 
only required to differentiate once with respect to time 
to find the gravitational wave strain. 

In figure [T3] we plot the cross-polarisation of the grav- 
itational wave strain as measured by an observer at 
d ~ 10 kpc. During the initial phase of the evolution, 
the strain is of the order of /ix ~ 10"^'', which repre- 
sents the lower limit of our numerical sensitivity. The 
strain amplitude is then excited during the nonlinear 
phase of the kink instability (i.e. after approximately 
50 — 60 ms) , at which point the strain approaches ap- 
proximately /ix ^ 10"^'*. 

As expected, the signal in figure [13] has the expected 
~ 1.8 kHz oscillations associated with the /-mode. In 
figure [14] we plot the Fourier transform of the cross- 
polarised strain, hx (/) for the entirety of the signal 
shown in figure [T3l (i.e. from t = — 380ms). The 
/-mode is clearly displayed here as a large peak, with 
a significantly smaller peak located at ~ 3.9 kHz. The 
signal amplitude present in a gravitational wave detec- 



tor is given by \/T hx (/) 



where T is the damping 



21, 



(16) 



time of the oscillation. The dominant mechanism for 
/-mode damping in neutron stars is through gravita- 
tional wave emission"^, which gives a damping timescale 
of T - 100 - 300ms [e.g. Ill, ^. For such a situa- 
tion, this gives a signal amplitude that is well below the 
detectable limit, even for third generation gravitational 
wave interferometers such as the proposed Einstein Tele- 
scope [78| . We discuss this in significantly more detail in 
section IIVI 

Figure [131 also exhibits large amplitude oscillations at 
significantly lower frequencies than the /-mode. These 
can be seen in the temporal evolution oi hx (figure I13p , 
and more clearly in the Fourier transform of this signal 
(figure [TJ]) . For our fiducial simulation presented here, 
the amplitude of the gravitational wave signal in these 
lower frequency modes is as large, if not larger, than the 
amplitude of the /-mode signal. The dominant damping 
mechanism for these modes is largely unknown, although 
it is generally expected that they will last significantly 
longer than the /-mode signal. Given that the signal 
amplitude scales as VT, these Alfven modes could be 
significantly more detectable than the /-modes given a 
large-scale magnetic field rearrangement following a mag- 
netar fiare. This is especially pertinent given that the fre- 



where an overdot denotes a time derivative, d is the dis- 
tance between the observer and the source and lij is the 
reduced quadrupole moment given by 

/ p i XiXj - -Sijx'^ j dV. 



i-ij — 



(17) 



Rather than numerically evaluating lij at each time-step 
and differentiating twice to get the gravitational wave 



^ Note that figure [T3l exhibits no damping of the /-mode through- 
out the long evolution time. This is because we are working in 
the Cowling approximation, which implies the system does not 
lose energy to gravitational radiation. A more thorough study 
involving full general relativistic effects is warranted, however 
we note the extremely large computational costs of simulations 
evolving the full spacetime render such a task presently intangi- 
ble. 
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FIG. 12: Equatorial slice of the absolute value of the magnetic field for four different time-points in our fiducial simulation - 
i.e. with average magnetic field Bis = 13, corresponding to an Alfven wave crossing time of 5ms. The magnetic field strength 
is plotted in a log-scale, with the maximum (shown in red) as |i3| = 4.8 x 10^^ G. Blue represents the minimum value of the 
magnetic field, and hence the neutral line is clearly visible in figures a, b and c. The thick black line is the surface of the star. 
The figures are a) i = ms, b) i = 25 ms, c) t — 50 ms and d) i = 195 ms. 



quency of these modes coincide with LIGO and VIRGOs 
most sensitive frequency band. 

It is possible that the low frequency Alfven modes seen 
in the present set of simulations could somehow be associ- 
ated with the observed quasi-periodic oscillations in the 
tails of giant flares. If this were true it would imply a 
damping time longer than 10s of seconds, which would 
significantly enhance the possibility of detection of grav- 
itational waves. Having said that, we have no direct ev- 
idence that the modes we see are in anyway associated 
with the QPOs, other than them existing in a similar fre- 
quency band. A more rigorous analysis of these modes is 
required to fully understand their relationship with ob- 
servations. However, an explicit analysis of these modes 
is proving elusive due to the violent dynamics associated 
with the magnetic field rearrangement. Moreover, it is 



expected that the spectrum of these mode will be ex- 
tremely dense (if not form a mathematical continuum), 
significantly hampering our efforts to provide a detailed 
expose on the existence of these modes. Therefore, for 
the moment, we shall have to be satisfied with the de- 
tailed literature analysing Alfven spectra of magnetars 
in the linear regime [e.g. dSt, [20, :22, :23fc ^ ^, i54i, i29|] . 



D. Quasi-Equilibria 

Figures [TUH and [T^Jl show a typical quasi- equilibrium 
snapshot a.t t — 195 ms — 39 ta- In some ways this snap- 
shot resembles a "twisted torus" configuration seen in 
the nonlinear evolutions of Braithwaite 149;] , Braithwaite 
and Spruit [ |8Q |, and in numerous semi-analytic equilib- 
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FIG. 13: Cross-polarisation of the gravitational wave strain, 
hx, from our fiducial simulation as measured by an observer 
at 10 kpc. The nonlinear phase of the kink instability induces 
a prominent gravitational wave signal which is present after 
approximately 50 ms. 
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FIG. 14: Fourier transform of the cross-polarisation of grav- 
itational wave signal, /ix(/), as a function of oscillation fre- 
quency for the fiducial simulation. The fundamental /-mode 
can be seen as the sharp peak at approximately 1.8 kHz. The 
modes present at less than 300 Hz are the lower frequency 
Alfven modes. The density of this part of the spectrum can 
be clearly seen in the plot inset, which shows a zoomed in 
region between and 500 Hz. 



Each of our figures showing time evolution of global 
quantities (i.e. figures [SJ IHl HH and [T^ show residual 
kinetic motion remaining in the star following an evolu- 
tion lasting almost 80 Alfven crossing times. We term 
the situation we have reached a pseudo-equilibrium as 
we can not formally prove whether this state is close to 
a stable equilibrium or not. It is clear that the global 
structure of our pseudo-equilibrium does not change fur- 
ther with time, in that all global quantities are tightly 
bounded. This is a necessary condition for claiming sta- 
bility, however it is not sufficient. Given the recent re- 
sult of Lander and Jones [50| , who has shown that vari- 
ous twisted-torus configurations originally expected to be 
stable [e.g.|43,[5H actually show instabilities, it is becom- 
ing increasingly likely that the gamut of stable magnetic 
field equilibria form an extremely complicated set. 

In the exterior region of our star we evolve the mag- 
netic field according to the same prescription as outlined 
in section |TT1 however we impose a low density atmo- 
sphere (patm ~ 10~^Pc where Patm and pc are the atmo- 
spheric and central density respectively). We note that 
this is not a good model for a magnetosphere as we do 
not impose a force-free condition. Moreover, we include 
boundary conditions that are close to the surface of the 
star and do not allow significant evolution of the mag- 
netic field (i.e. we use Dirichlet boundary conditions at 
the outer bounding box of our simulations). However, we 
do evolve the magnetic field in the stellar exterior, which 
allows us to determine the extent to which the external 
field is altered by the magnetic field instability inside the 
star. 

Figure [T51 is equivalent to figure [TUH - however showing 
the full domain of the Cartesian grid and not truncating 
the magnetic field lines at the surface of the star. Despite 
the large-scale rearrangement of the internal magnetic 
field topology, one can see from this figure that the exte- 
rior magnetic field remains largely unchanged throughout 
the nonlinear saturation of the instability. This can be 
clearly seen by remembering that the initial field is purely 
poloidal, which is in close approximation to the exterior 
part of the field shown in figure 1151 



E. Effect of Magnetic Field Strength 



rium derivations that include both poloidal and toroidal 
field components. Indeed, the left hand side of the star 
as seen in figure [TUH is well approximated by a twisted- 
torus (i.e. with toroidal components of the field confined 
to the closed field lines of the poloidal field), and the inte- 
rior is threaded by a dominantly poloidal field. However, 
the remainder of the star exhibits large non-axisymmetric 
structures with toroidal and poloidal field apparently in 
equal abundance. This non-axisymmetric structure can 
clearly be seen in figure [T2H . which also highlights the 
fact that the neutral line no longer lies on the equatorial 
plane, and has been severely disturbed by the reconfig- 
uring of the field. 



Until now we have provided a somewhat comprehen- 
sive analysis of a single, fiducial simulation with po- 
lar surface magnetic field strength of Bi^ — 8.8, cor- 
responding to an Alfven crossing time of ta — 5.0 ms. 
Herein we present the same simulation with differing ini- 
tial magnetic field strengths between polar surface fields 
of i3i5 — 6 and Bi^ = 55, corresponding to Alfven cross- 
ing times of ta — 7.3 ms and ta — 0.9 ms respectively. 
It is worth noting that our lowest field strength simula- 
tion is now only a factor of about two stronger than that 
observed in magnetars. 

In figure \T^ we plot the evolution of the m = 1 com- 
ponent of Cm (Be/,) defined by equation (|15p . We have 
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FIG. 15: Same as figure llOt l. however now showing the full 
domain of the grid, and not truncating field lines at the sur- 
face of the star, which is shown as the green circle. Despite 
a significant rearrangement of the magnetic field inside the 
star, the exterior has remained significantly unchanged, as 
evidenced by the fact that it is almost purely poloidal. 



normalised the temporal axis now to the Alfven crossing 
time of the system. One can trivially see an invariance 
in the growth time of the instability as a function of the 
Alfven crossing time, and hence magnetic field strength. 
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FIG. 16: (Top panel) Evolution of Cm=i (B0) as defined in 
equation (|15p as a function of time for our fiducial F = 2 
model with different magnetic field strengths. The tempo- 
ral unit is normalised to the characteristic Alfven timescale 
throughout the star to show invariance of the instability 
timescale. (Bottom panel) Poloidal-to-total magnetic energy, 
Ep/E, as a function of time for our fiducial F = 2 model with 
different magnetic field strengths. Despite remaining residual 
motions in the system after up to 100 Alfven timescales, all 
models quasi-equilibrium state have Ep/E ~ 0.65. 



IV. EFFECT OF EQUATION OF STATE 



In the bottom panel of figure[T6]we plot the poloidal-to- 
total magnetic energy ratio as a function of time, where 
we have again normalised the temporal unit to the char- 
acteristic Alfven crossing time of the system. While the 
instability growth timescale of the system can not be 
as readily established from this plot, one can still see 
here that this timescale trivially scales with the Alfven 
crossing time of the system. Moreover, we see that the 
poloidal-to-total magnetic energy ratio of the equilibrium 
state is independent of the magnetic field strength of 
the simulation. That is, regardless of the magnetic field 
strength, the final quasi-equilibrium has Ep/E ~ 0.65. 



In Zink et al. (32| we established a relationship between 
the surface magnetic field strength and the gravitational 
wave emission (both energy and strain). In the present 
article we generalise this result to include different poly- 
tropic equations of state, and hence stellar radius and 
mass dependencies. Rather than review those previous 
results here, we give a complete analysis in the following 
section. 



The primary purpose of this section is to generalise 
our gravitational wave results of Zink et al. 32] to in- 
clude the dimensions of the star. That is, in Zink et al. 
[32| we evaluated the gravitational wave emissions due to 
a reconfiguring magnetic field, induced by the magnetic 
field instability. This is deemed to be a toy model for the 
aftermath of a magnetar flare, in which the flare causes 
some form of rearrangement of the magnetic field. Note 
that we are not saying that the kink instability causes 
a magnetar flare, rather that we are simply determin- 
ing the energy conversion between magnetic field, fluid 
and spacetime dynamics utilising the kink instability as 
a mechanism to generate such motions. 

In Zink et al. |32| we determined a power-law relation- 
ship such that the gravitational wave strain is propor- 
tional to the surface magnetic field strength to the power 
3.3. This only took into account our fiducial F — 2.0 
polytrope, and therefore was not determined as a func- 
tion of any other stellar parameters. In this section we 
repeat the analysis using a range of soft and hard EoSs 
to determine a scaling law relating the gravitational wave 



14 



Model 


r 


K 


M 


R 


Freq. 








[Mo] 


[km] 


[kHz] 


AO 


1.67 


0.0400 


1.31 


19.28 


1.2 


BO 


2.00 


0.0269 


1.31 


12.68 


1.7 


CO 


2.34 


0.0195 


1.31 


11.54 


1.8 


DO 


2.46 


0.00936 


1.31 


8.47 


2.5 



strain, surface raagnetic field strength, stellar radius and 
mass. 



TABLE I: Equations of State parameters used throughout the 
article. Note that model BO is our non-rotating fiducial model 
from [3ll.l3^1 and section Hill and models CO and DO are EoSs 
II and A respectively. The label represents no rotation, M 
is the gravitational mass, R the equatorial radius and Freq. 
in the final column is the fundamental /-mode frequency of 
the system giving the characteristic fluid timescale. 




FIG. 17: Maximal strain, /i™'"', as a function of polar surface 
magnetic field strength for non-rotating models presented in 
table m These models all have constant stellar radius. 



Table U shows the EoS parameters that we utilise 
herein. For each EoS we have constructed a series of 
models with central magnetic field 1.6 x 10^^ G < |i?c| < 
2.7 X 10^^ G. Depending on the EoS, these models will 
have different surface magnetic fields and different Alfven 
timescales. For each model we plot the maximum value of 
the cross-polarisation of gravitational wave strain, h'^^^, 
as a function of the surface magnetic field strength at the 
pole, i?poie, in figure flTl As the star becomes more com- 
pact (i.e. as r increases), the strain evaluation for the 
models with weaker magnetic fields becomes more diffi- 
cult. Therefore, for the below analysis, we restrict our 
attention only to the stronger field strength models for 
which the power-law relation is valid. 

Figure [T7] shows the dependence of strain on the 
mangetic field and also the equation of state, where ex- 
plicit dependence is imposed on the radius of the star. 
The full functional dependence of this relation will also 
include the stellar mass (or alternatively the central den- 
sity). We have therefore produced extra series' of models 
for equations of state BO and CO, allowing the mass to 
vary. This allows us to perform a full, least squares fit 
over the four parameter space h'^'^^, R, M and Spoio, 
finding the following relation 
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(18) 



This result is consistent with that presented in Zink et al. [32] in that ft,™'"' goes approximately to the third power 
of the magnetic field strength. The additional factors in this equation are the dependence on the stellar radius and 
mass, R and Af , which scale almost to the fifth and second powers respectively. 

As previously discussed, a majority of the energy in the signal associated with ft™^'' is in the f-mode. We can 
therefore calculate an approximate amount of energy emitted in gravitational radiation by noting |81| 



27r2d2/2c3 
^GW ^ / {h") dt, 



(19) 



where / is the /-mode frequency of radiation. Assuming a gravitational wave damping time of 100 ms, we find a 
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power-law relation for the total energy emitted in gravitational radiation via the /-mode to be 
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(20) 



r 



This highly non-linear relationship with the radius 
could have implications for gravitational wave detection. 
If nature is kind to us, and we find that the EoS of neu- 
tron stars are soft, implying their radii are large with 
respect to the fiducial model, then gravitational wave de- 
tection of /-modes from magnetar flares becomes more 
likely. 

It is worthy of note that, while the lower field values 
for the DO model are off-kilter, and hence have been ex- 
cluded from the above analysis, the three stronger field 
models are consistent with the above relation. We have 
confirmed this by removing all of the DO models, at which 
point the least squares analysis gives similar results to 
those presented above. 

In the above we have calculated the maximal strain 
seen in our simulations rather than a time integrated 
signal such as a traditional root-sum-square amplitude 
^rss = S^oo (^^) ^^ [^-S- ISa|- The main reason for this 
is that a time integrated signal would require us to de- 
fine a time for the onset of the instability. This is be- 
cause there is no intrinsic damping in the system (i.e. we 
are working in the Cowling approximation, implying the 
dominant /-mode damping mechanism is not present), 
implying we would have to ensure we only integrate for 
T seconds (where T is the damping time of any particular 
mode) following the onset of the instability. Integrating 
beyond this time would artificially grow the hj-ss signal. 
Moreover, were we to integrate from the very beginning 
of our simulation out to some fixed time, this would in- 
troduce a systematic error as a function of the magnetic 
field strength, due to the onset of the instability being a 
function of the Alfven timescale of the system. 

Instead of the aforementioned approach, we calcu- 



hif) 



for our simula- 



late the signal amplitude, 

tions as a function of the frequency, such that this 
value can be directly compared with the noise power 
spectral density, y/\Sh (/)|, of individual gravitational 
wave detectors, and the amplitude signal-to-noise ratio, 

y/T h{f) /\/\Sh (/)|, can simply be read of the resul- 
tant figure. To calculate our signal amplitude, we take a 
Fourier transform of a portion of our derived, hx that is 
post kink instability saturation and lasts approximately 
150 ms (this value ensures we get a reasonable number 
of oscillations in the lower portion of the spectrum) . We 
then multiply by the relevant damping times and plot 
the results in figure fTSl 

In particular, each panel of figure [18] represents dif- 

* We note the error in our published article Zink et al. |33l where, 
in equation (2), we have an anomalous dependence on the dis- 
tance to the source in the relationship for the gravitational wave 
energy. 



ferent EoSs given in table HI For each EoS we have run 
multiple models with different magnetic field strengths, 
and for each simulation have located the /-mode signal 
amplitude and the maximal Alfven mode signal ampli- 
tude. The /-modes (coloured boxes) are plotted in figure 
[T8]assuming a damping time between T = 50 and 200 ms. 
One can see that the /-mode frequency for the different 
EoSs ranges from ^ 1.2 kHz for models AO to ^ 2.5kHz 
for models DO. Moreover, for all of the models presented 
herein, only the largest strength magnetic field model 
(with surface field strength Bis = 18) with the softest 
EoS is observable by the Advanced LIGO detector. It is 
worth noting that the strongest magnetic field observed 
to date in a magnetar is almost a full order of magnitude 
less than that modelled for this particular star. The de- 
tectability situation of /-modes could change slightly if, 
for example, neutron stars have interior toroidal mag- 
netic fields an order of magnitude larger than the ob- 
served dipole poloidal field. Although this situation is 
unlikely due to stability arguments, a point we discuss in 
more detail below. 

For each EoS in figure [18] we also plot two or three 
full spectra with assumed damping time of T = 100 ms. 
While the /-mode is abundantly clear in all of these 
curves (except from model AO with B15 = 18 - a point 
we discuss below) , lower frequency Alfven oscillations are 
also clear in each of these simulations. As discussed pre- 
viously, the damping time for these modes is largely un- 
clear, but is likely significantly longer than a millisec- 
ond. Therefore, for each maximal Alfven mode, we have 
also plotted a vertical dumbbell (i.e. a skinny line with 
bulbous ends) with damping times between T = 10 ms 
and 1 s. For significantly lower values of magnetic field 
strength, these bars become closer to the lower limit of 
detection for both the Einstein Telescope and possibly 
even Advanced LIGO. This is especially true when one 
considers that the damping time of these modes could, 
in fact, be minutes or even significantly longer. These 
modes, as with the /-modes discussed above, are also 
better excited for softer EoSs, corresponding to larger 
neutron stars. 

It is important here to mention the signal amplitude 
of the AO model with B15 = 18. This model has 
an Alfven crossing time of approximately ta = 1.1ms, 
which corresponds to a fundamental Alfven frequency of 
/a — 910 Hz. Moreover, the /-mode frequency of the AO 
simulations is approximately 1170 Hz. One can therefore 
see from the Fourier transform that the region around 
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FIG. 18: Gravitational wave signal amplitude, \/r Ki (/) , against oscillation frequency for the four equations of state used 

herein - AO (top left panel), BO (top right), CO (bottom left) and DO (bottom right). For each box model there are several 
magnetic field strength simulations shown. The /-mode signal amplitude is represented by a coloured box at the appropriate 
frequency (shown in table |l]), where a source lasting between 50 and 200 ms has been assumed. The coloured vertical dumbbells 
(i.e. skinny lines with bulbous ends) represent the maximum Alfven mode seen in the Fourier transform. These assume a 
constant source lasting between 10 ms and 1 s. Both modes scale with vT, implying one can easily extrapolate to alternative 
values of the damping time. We have further plotted the entire spectrum (assuming a damping time of T = 100 ms for numerous 
models to give an indication of the spectra present. Finally, the noise power spectral density curves for LIGO, AdvLIGO and 
ET are shown in each panel. These have been taken from the review article of Sathyaprakash and Schutz JT^] . 



the /-mode frequency is not a clean peak, but shows a 
rather dense spectrum consistent with expectations from 
the lower field simulations. In fact, while there exists a 
distinct peak at the /-mode frequency, this is not the 
largest peak present in the Fourier transform. Rather, 
the largest peak, less than a factor two larger than the 
/—mode peak, is located at approximately 1.35 kHz. We 
have indicated this on the top left panel of figure [T8l 

Our conclusion from figure |T5] is the following: Assum- 
ing that a giant magnetar flare is somehow related to a 
catastrophic rearrangement of the core magnetic field, the 
gravitational wave signal associated with f -modes are not 
observable with present or near-future gravitational wave 
observatories. Lower frequency Alfven modes do, how- 
ever, provide an enticing alternative where efforts could 
be concentrated. Of course, the biggest hurdle to under- 
standing the gravitational wave emission associated with 



these Alfven modes is understanding their damping time. 
If they are sufficiently long-lived, then they may be de- 
tectable in the relatively near future. 



V. ROTATION 

Frieman and Rotenberg [S^] were the first to show that 
rigid-body rotation only has a significant effect on hydro- 
magnetic equilibria when the fluid-flow velocity is of the 
same order as the Alfven velocity. Certainly, in newly- 
born neutron stars this is expected to be the case in large 
regions of the star, even for magnetar field strength stars. 
Pitts and Tayler [82| studied rotating, toroidal magnetic 
fields in cylindrical geometries, showing that the pertur- 
bations can be countered by sufficiently large rotational 
velocities, implying the suppression of certain modes of 



17 



Model 


Be 


n 


flp/ae 


rn 


TA 




[10'" G] [Hz] 




[ms] 


[ms] 


BO 


100 





0.99 


oo 


2.4 


BlOO 


100 


100 


0.99 


10.0 


2.4 


B200 


100 


200 


0.97 


5.0 


2.4 


B300 


100 


300 


0.94 


3.3 


2.4 


B400 


100 


400 


0.90 


2.5 


2.3 



[T9l where we plot measures of the rotational energy in 
the system for five rotating models, all with initial spin 
period of to = 5.0 ms. In the top panel of figure [TOl we 
plot the angular momentum, J, defined as 



(21) 



TABLE 11: Rotating models. Our model has equation of state 
B given in table HI with central field of Be = 1.0 x 10^'' G. 
Moreover, Q is the rotational frequency, ap/ue is the ratio of 
equatorial to polar radii and rn & ta are the rotational period 
and Alfven timescale respectively. 



instabilities. 

The first work on the stability of purely poloidal, ro- 
tating fields was from Geppert and Rheinhardt [53], who 
performed three-dimensional, non-linear numerical simu- 
lations by use of a spectral code. They found that stars 
rotating with sufficiently high rotation speeds, quantified 
as Ha/^ ^0.1, and with roughly aligned magnetic and 
rotation axes, will have suppressed instabilities. In con- 
trast, Braithwaite [83] also performed nonlinear, global 
simulations of rotating poloidal magnetic fields, finding 
that the initial linear phase of instability growth was 
not affected by the presence of rotation. However, while 
Braithwaite |83| found that the nonlinear phase is af- 
fected by rotation, the instability was always present re- 
gardless of the rotation speed. Lander and Jones [5J] re- 
cently performed linear simulations of rotating poloidal 
fields, showing that some, but not all modes were sta- 
bilised by the presence of rotation. 

Herein we attempt to resolve some of the afore- 
mentioned contradictions, simultaneously performing the 
first simulations of rotating poloidal fields in general rel- 
ativity. 

For our rotating simulations we revert back to the fidu- 
cial EoS - i.e. a r = 2 polytrope as described in section 
mil We have created a series of rotating models utilising 
the spectral solver lorene with central magnetic field 
strength of Be = 1.0 x 10^^ G, corresponding in the non- 
rotating limit to a surface field strength of i?i5 = 16. We 
summarise the properties of this models in table |TT1 



A. Effect of Magnetic Field on Rotation Rate 

Our simulation method is identical to the previous sec- 
tion. Initial conditions are created on a spectral grid us- 
ing LORENE, which is then mapped to our Cartesian grid 
measuring 120 grid-points in all directions. No perturba- 
tion is added, and we allow the simulation to evolve. 

It is worth noting from the outset that we see a not- 
so insignificant spin-down of the rotating star due to the 
presence of the magnetic field. This is shown in figure 



J= / UaT^^d^X, 



where n^ is the four- vector that is hypersurface orthogo- 
nal to surfaces of constant time, 7ij is the three-metric of 
the spatial hypersurface E and T'"^ is the stress-energy 
tensor. In the bottom panel of figure[Tn]we plot the ratio 
of the rotational kinetic energy to gravitational binding 
energy T / \ W\- For definitions of these quantities see 
Stergioulas [83]. It is worth noting that, for our calcu- 
lation of the angular momentum, we assume the system 
remains axisymmetric. i.e. that d^ is a Killing vector. 
While this will introduce an error into the calculation 
of the angular momentum, we expect this to be minimal 
due to the fiuid remaining almost axisymmetric, with de- 
viations from this only arising due to non-axisymmetries 
induced by the magnetic field. 

Model A in figure [T51 has zero magnetic field through- 
out the star. The angular momentum and T/]VF] are 
extremely well conserved over the twenty spin periods 
shown in this figure. We have evolved such a simula- 
tion for almost 700 ms, showing conservation of angular 
momentum on the order of 7%. 

Models B - E of figure [H are all B200 from table HH 
model B has 120^^ grid points with the outer boundary 
located at approximately 1.4 times the surface of the star 
(i.e. this is our canonical simulation set-up that will be 
used throughout the remainder of the article) , model C 
has 180"^ grid-points, however the outer boundary of the 
star has been moved to a greater radius such that the 
grid resolution across the star is the same as in model 
B. In model D we have 150^ grid-points with all other 
quantities being the same as model B. Finally, in model 
E we have used 150'^ grid-points and implemented a linear 
extrapolation boundary condition for the magnetic field 
at the outer-boundary of the domain. This is in contrast 
to the Dirichlet boundary conditions used for all other 
simulations. 

One can see from figure [19] that the loss of rotational 
energy is almost independent of our simulation method 
(in terms of resolution, boundary location and magnetic 
field boundary conditions). In fact, the dominant factor 
in the loss of angular momentum is the strength of the 
magnetic field, implying this scales almost linearly with 
the Alfven timescale of the system. We note here that 
model E in figure [TOj has a slight up-turn in angular mo- 
mentum and T/]VF] after approximately 85 ms. This is 
not a conservation of rotational energy, rather this sim- 
ulation develops a numerical instability that develops at 
the boundary of our domain and causes the simulation 
to crash shortly after 100 ms. Higher-order interpolation 
schemes at the boundary could act to stabilise such a 
scheme, however figure [TOj indicates that this would not 
change the rate of rotational energy loss from the system. 
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FIG. 20: Normalised absolute value of the angular momen- 
tum, I J|/|Jo|, where Jo = J{t = 0), as a function of time for 
our models in table HH 



FIG. 19: Angular momentum, J, and r/|iy| for models with 
initial rotational period of 200 Hz. Model A has 120^ with zero 
magnetic field. Models B - E are all B200; model B has 120^ 
grid-points (i.e. is model our canonical simulation set-up), 
model C has 180'^ grid-points with outer boundaries located 
further from the star (such that the resolution across the star 
remains the same as model B). Model D has 150^ grid-points 
with the outer boundary located the same distance as in B, 
and E has 150^ grid-points with linear extrapolation for the 
magnetic field used at the outer boundaries rather than the 
standard Dirichlet boundary conditions for the other models. 



Because of the loss of angular niomentum in our system 
we are presently restricted to only discussing the onset 
and nature of instabilities, and can not provide any in- 
sight into the nature of pseudo-equilibria attained after 
long evolution times. In figure [20] we plot the absolute 
value of the angular momentum, normalised to the initial 
state, as a function of time for our models presented in ta- 
ble ini One can see here that our simulations retain more 
than 50% of their angular momentum for approximately 
100 ms, corresponding to more than 40 Alfven timescales 
and between 10 and 40 initial rotational periods. As we 
shall see below, these timescales are long enough to dis- 
cuss much of the dynamics of the systems in terms of the 
evolution of the varicose and kink instabilities. However, 
we stress that these losses of angular momentum imply 
that the end states of our simulations are not necessarily 
representative of equilibrium states of rotating neutron 
stars. We are therefore reticent to discuss such equilibria 
in the present article, keeping our discussion to that of 
the dynamics of the instability. 



B. Instability 

In figurel^we plot Ci (B^) and Ep/E (top and bottom 
panels respectively) for the series of models presented in 



table mi with initial rotational frequencies of between 
and 400 Hz. 

If we focus our attention on the initial first few ms of 
the simulation, particularly in terms of the Ep/E quan- 
tity (lower panel), we see a large trough reaching as low 
as Ep/E ^ 0.8 for the 400 Hz model, but smaller for lower 
rotation rates. This is a result of the mapping between 
the spectral grid and our Cartesian grid. The effect is to 
introduce strongly toroidal components into the magnetic 
field in the first few milliseconds, particularly around the 
equatorial region near the neutral line of the field. As can 
also be seen in this figure, this non-physical effect also 
vanishes after the first few Alfven crossing timescales, 
and the simulation reduces to an almost purely poloidal 
state. 

The immediately striking part about figure [21] is that 
the instability timescale in terms of the Cm {B^) quan- 
tities is independent of the rotational velocity. For all 
initial rotational velocities, we see the Ci (B^) grow ex- 
ponentially by over six orders of magnitude in approxi- 
mately 50 ms. This is the same growth timescale of the 
instability in terms of the poloidal to total energy ratio 
for the case with zero rotation rate. However, this is not 
the case for the rotating simulations. For example, if we 
focus our attention on model BlOO (blue dashed line), 
one can see that the timescale for the ratio of poloidal 
to total magnetic energy to reach the canonical "equi- 
librium" value of Ep/E ~ 0.65 is approximately 80 ms. 
In this case, the rotation has slowed the development of 
the instability in terms of magnetic field reconstructions. 
This is more extreme for the 300 and 400 Hz models (dot- 
ted green line and dot-dash pink line respectively) , which 
are shown to be stable for the full 115 ms evolution shown 
here. 

The seeming contradiction in growth timescales can be 
understood by noting that the CmB^ quantities trace any 
change in the structure of the magnetic field around the 
equator, while Ep/E only tracks the growth of toroidal 
magnetic field. In terms of the two known instabilities 
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FIG. 21: Magnetic field instability in terms of both Ci (B^) 
(top panel) and poloidal-to-total magnetic energy ratio, 
Ep/E, for the models presented in table [Til The individual 
lines are the same as those for figure [20l with the black line 
representing the non-rotating model. 



expected to act on such magnetic field configurations, 
CmBfj, is sensitive to both the varicose and kink modes, 
while the poloidal to total magnetic energy ratio is only 
sensitive to the kink mode. This can be understood as 
the varicose mode only acts to change the cross-sectional 
areas of flux tubes around the equator, keeping them or- 
thogonal to equipotential surfaces. The kink instability, 
on the other hand, acts to directly increase the amount 
of toroidal magnetic field in the star. 

To highlight the aforementioned affect we look at 
three-dimensional plots of the magnetic field lines to ex- 
pose the global dynamics of the simulation. In particular, 
in figures !^ we show three-dimensional snapshots of the 
BlOO simulation (i.e. corresponding to the dashed blue 
line in figure I^Tj) . The four snapshots in this figure are 
shown after times t — ms (figure [22k). t — 17 ms (figure 
[22b), i = 27 ms (figure [22t), i = 42 ms (figure [22H). 

Figure [52k represents a similar initial condition to that 
seen in figure [TOk . In figures [22b and [22b we see the de- 
velopment of a strong varicose mode that has disrupted 
the magnetic field in a significantly more catastrophic 
manner than that seen for the non-rotating simulations. 
This is especially seen when these figures are compared 
with figure [TUb . If we compare the time of the three- 
dimensional snapshots seen in figures Wlb and [22b with 
the instability plots given in figure [21] we see that these 
represent points at which the Cm (B^) quantity is grow- 
ing exponentially, while Ep/E remains almost constant. 

As discussed in section [Til Bl the varicose mode was ex- 



cited most strongly in the non-rotating case as an m = 4 
mode. Although the presence of the varicose mode is 
a physical effect, the dominance of the m — A mode is 
attributed to excitations caused by the Cartesian grid. 
The present rotating cases however, are moving with re- 
spect to the stationary grid, implying this mode is not 
preferentially excited. Moreover, we see these modes ap- 
proximately equally in all azimuthal wavenumbers. We 
therefore conclude that the varicose mode seen in our 
simulations are completely physical. Moreover, from the 
top panel of figure[2Tl we conclude that the varicose mode 
is, in fact, unstable. 

Finally, figure [22b l. shows a strong kink instability 
developing, corresponding to a transference of poloidal 
magnetic field energy to toroidal. As mentioned, this 
figure is shown after i = 42 ms which, in figure [21] corre- 
sponds to Ep/E ^ 0.85. 



In figures [23| we again plot three-dimensional snap- 
shots, however this time for the rotating B400 model. 
The snapshots shown here are a.) t ~ ms, b) i = 27 ms, 
c) i = 73 ms and d) i = 148 ms. In figure [23k we see 
a slight change to the initial conditions to that seen in 
figure[Tni In particular, the red magnetic field lines cover 
a larger area. This is due to the initial rearrangement of 
the magnetic field, which also acts to push the neutral 
line to a larger equatorial radii. As the field lines are 
seeded at a constant radii, the visual effect is that these 
field lines cover a larger fiux surface. It is worth noting 
that these field lines are still wholly contained within the 
star. 

Once again, figure [23b exhibits a strong varicose mode 
growing early in the simulation. This motion remains 
unchanged for approximately 50 more ms. After ap- 
proximately 70 ms one sees Ep/E ~ 0.95. The corre- 
sponding snapshot for this figure exhibits a single kink 
around the star, which rotates with the star, remaining 
virtually unchanged for approximately another 100 ms. 
It is extremely tempting to say that this is an equilib- 
rium configuration for a fast rotating (i.e. tq ~ r^), 
strongly magnetised star, as figure [23H. shown after al- 
most 150 ms, exhibits a very similar structure. Almost 
180 ms into the simulation one finally sees the kink in- 
stability go unstable, and Ep/E evolves to the canonical 
model of Ep/E ~ 0.65. We have not shown this in figure 
[2T]as we believe it could possibly be mis-leading due to 
the substantial rotational energy loss in the star. It is 
therefore not clear from this investigation whether the 
kink instability is only present in these simulations be- 
cause enough angular momentum has been lost from the 
simulation such that the condition ta > tq is once again 
satisfied, and the kink instability can act. 



It is worth pointing out that we have only explored 
models for which ta ^ tq. We have not yet broached 
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FIG. 22: Time evolution of rotating model BlOO defined in table|lT]- i.e. with central magnetic field Be = 1.0 x 10^^ G, initially 
rotating at 100 Hz. This gives tci/ta ~ 4.2. The figures are a) t — ms, b) i = 17 ms, c) f = 27 ms and d) t — 42 ms. A 
description of the figure is given in the caption to figure ITOl 



models whereby r^ » ta, implying the rotational ve- 
locity is dominating over the Alfven timescale. It is in 
this regime where one expects from linear analysis and 
Newtonian simulations [53| that the kink instability will 
be completely suppressed. 



VI. CONCLUSION 

We have performed three-dimensional, general rela- 
tivistic ideal MHD simulations of rotating and non- 
rotating polytropic neutron stars with initially purely 
poloidal magnetic field geometries. This was accom- 
plished by evolving initially self-consistent solutions of 
the Einstein-Maxwell field equations under the ideal 
MHD assumption and the Cowling approximation. Par- 
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FIG. 23: Time evolution of rotating model B400 defined in table|lT]- i.e. with central magnetic field Be = 1.0 x 10^^ G, initially 
rotating at 400 Hz. This gives tq/ta ~ 1.1. The figures are a.) t — 0ms, h) t = 27ms, c) t — 73ms and d) t — 148ms. A 
description of the figure is given in the caption to figure ITOl Note that the field lines are seeded in the same location as figure [TOl 
however the fast rotation rate pushes the neutral line further from the center of the star, explaining the larger region covered 
by the red field lines. 



ticularly in the non-rotating case, purely poloidal fields 
are intrinsically unstable to the kink instability, which 
leeds to a catastrophic rearrangement of the magnetic 
field. Such reconfigurations allowed us to explore three 
separate phenomena: 

i) the nature of the instability itself, 



ii) the equilibrium configurations derived as steady- 
state solutions and 



iii) the gravitational wave emission from such magnetic 
field rearrangements which we interpret as a phe- 
nomenological model describing gravitational wave 
emissions from magnetar flares. 
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A. Hydromagnetic Instabilities 

The first part of this work associated with the instabil- 
ity and the pseudo-equihbria configurations represents a 
significantly more detailed presentation of the first non- 
linear GRMHD simulations of such instabilities 311. In 



particular, we showed that the kink instability acts on 
timescales associated with the Alfven crossing time of the 
star, TA, as expected from linear and non- linear analyses 
alike, saturating after approximately ISita- The kink in- 
stability acts on gravitational equipotential surfaces, and 
dominates near the neutral line of the poloidal field. This 
acts to induce strong toroidal components of the mag- 
netic field, which eventually saturates with a poloidal-to- 
total magnetic energy ratio of Ep/E ^ 0.65. 

We have made a first attempt at including rotation into 
our simulations, making these the first rotating GRMHD 
simulations with non-zero magnetic field in the atmo- 
sphere of the star. This provides a mild technological 
challenge due to the significant numerical difficulties as- 
sociated with both the stellar boundary and the exterior 
atmosphere region, particularly with reference to the re- 
covery of primitive variables transformations. 

There has been debate in the literature as to the affect 
rotation has on the varicose and kink instability. Gep- 
pert and Rheinhardt [53| found that sufficient rotation 
can act to suppress such instabilities, while Braithwaite 
[83| showed that rotation does not affect the presence or 
linear growth rate of such instabilities, although he only 
performed simulations of rotating toroidal fields. Our 
rotating simulations showed that rotation acted to sep- 
arate the varicose mode from the kink mode, with the 
former having a growth timescale unaffected by the pres- 
ence of rotation (i.e. still of order the Alfven crossing 
time) , while the latter being slowed by the rotation. De- 
spite the slowing of the kink mode, this mode was still 
present in all of our rotating simulations. It is worthy of 
note that our rotating simulations only begin to broach 
the regime in which tq ^ ta, but does not yet extend 
into the regime in which tq >> ta- 



torii, in that they contain predominantly poloidal re- 
gions, where the closed field lines are threaded by toroidal 
components of the field - for example figures[T0h.[22H and 
123d . In some sense, our equilibria models are roughly 
a cross between the non-axisymmetric configurations of 
Braithwaite [5i|, and the same authors axisymmetric 
equilibrium configurations [49j. The difference between 
these two equilibrium derivations was predominantly as- 
sociated with the position of the neutral line in the star. 

It is worth reiterating that the equilibria derived from 
our evolutions still show considerable kinetic energy in 
the system, even after up to a hundred Alfven crossing 
times. We describe this situation as a psewrfo-equilibrium 
as we can not formally prove how close this situation is 
to a stable, stationary equilibria. As discussed in the 
text, it is abundantly clear that the global structure of 
our pseudo-equilibria does not change further in time, 
and all globally measured quantities are tightly bounded. 
For example, the ratio of magnetic energy in the poloidal 
component of the magnetic field to the total magnetic 
field is Ep/E - 0.65. 

The stability of various twisted-torus configurations 
has recently been called into question [50]. In that 
paper, axisymmetric solutions of the baratropic Grad- 
Shafranov equation with mixed poloidal and toroidal 
magnetic fields were linearly perturbed, with all result- 
ing e quili bria found to be unstable. The Lander and 
Jones [50| 'type-I' solutions are those resembling twisted- 
torus models in that they have an interior toroidal field 
threading only the closed field lines of the poloidal field. 
Those authors however, could only probe equiblira so- 
lutions with Ep/E > 0.945^. The stability of equilib- 
rium configurations such as those presented herein and 
in Braithwaite [43, [51| is therefore still an open question 
- although, given the bounded nature of any kinetic mo- 
tion, we expect these pseudo-equilibria to be stable to 
arbitrary perturbations. 



C. Gravitational Waves from Stellar Oscillations 



B. Equilibria 

A largely unresolved field of research is that of de- 
termining the set of stable MHD equilibria possible in 
barotropic stars. Our evolutions push further the bound- 
ary of understanding. In particular, our initially ax- 
isymmetric systems invariably evolve non-axisymmetric 
structures early into their evolutions. All of our subse- 
quent equilibria are therefore likely to retain such non- 
axisymmetric structures, implying we are not likely to 
converge on any of the semi-analytic equilibrium mod- 
els used in the literature [e.g. [a, |43, [S^, nor are we 
likely to converge on the axisymmetric numerical results 
of Braithwaite and Nordlund [4^ , Braithwaite [431 ■ ^^~ 
spite this, many of the end-states that we derive contain 
considerable portions of the star that resemble twisted- 



A further application of the simulations presented 
herein is to understand the gravitational wave emission 
from magnetar fiares. We reiterate that we are not advo- 
cating that kink instabilities are responsible for magnetar 
fiares. Instead, we are utilising the kink instability to 
mimic a global reconstruction of the internal magnetic 
field expected immediately following a magnetar fiare. 
This is pertinent given the current experimental searches 



^ The quoted maximum value in Lander and Jones |50| is Ex/E = 
0.045, where Ex = E — Ep is the magnetic energy in the toroidal 
component, however this has been calculated as a volume integral 
that includes the exterior region of the star. The same calculation 
restricted to the region bounded by the star gives Ex/E = 0.055 

[8i. 
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for gravitational waves from magnetar flares from the 
LIGO Scientific Collaboration [87,-M]- 

Our key results in the process of determining overall 
gravitational wave emission are equations ^TE\\ and (PUI . 
In particular, equation (llSp shows the expected gravi- 
tational wave strain amplitude given a star of mass M, 
radius R, polar magnetic field strength -Bpoic and located 
dkpc from Earth. In Zink et al. 32] we presented a first 
self-consistent estimate of this given a single equation 
of state, finding that the gravitational wave strain scales 
with approximately the third power of the magnetic field. 
The present paper extends this to show that the gravi- 
tational wave strain scales with almost the fifth power of 
the radius and the mass squared. 

To derive the previous result we explored multiple 
polytropic equations of state, ranging from polytropic 
index of F = 1.67 to F = 2.46. The conclusion for all 
of these simulations was the same - the gravitational 
wave signal associated with f -modes is not observable 
with present or near-future gravitational wave observa- 
tories, including the proposed Einstein Telescope. Of 
course this result comes with the caveat that we are not 
observing extremely bloated stars with larger than ex- 
pected radii. 

Another factor that could alter our conclusions regard- 
ing the /-mode excitation is that we have only explored 
motions due to unstable poloidal fields. As discussed 
above, various mixtures of poloidal and toroidal field con- 
figurations could exist in the interior of neutron stars, 
including configurations with a relatively strong toroidal 
component. Such configurations would have stronger 
magnetic fields, hence higher coupling with the /-mode 
and consequently greater gravitational wave luminosities. 
Whilst an exploration of these effects is certainly war- 
ranted, we do note that such fields are unlikely to pro- 
vide the orders-of-magnitude difference required in grav- 
itational wave emission to become detectable. 

Although the previous /-mode result represents a 
somewhat disappointing science case for gravitational 
wave observatories, there is potential for excitement in a 
somewhat unexpected region of the spectrum. Our sim- 
ulations have shown the presence of strong oscillations 
in the low, i.e. lO's to lOO's of Hz, region of the spec- 



trum. These oscillations are consistent with the propa- 
gation of Alfven waves in the interior of the star. Our 
results only show the indication of these modes, how- 
ever the detectability requires a better knowledge of the 
various damping mechanisms relevant for these modes 
of oscillation. In figures [TH] we have shown these modes 
assuming a damping time of 10 nis to 1 s. The gravita- 
tional wave signal amplitude, and hence detector sensi- 
tivit y, s cales as vT, where T is the length of the signal 
[e.g. |78|. Therefore, an increase in damping time from 
1 s to 1 hr increases the signal amplitude by a factor of 
60. To the best of our knowledge, estimates of expected 
damping times for Alfven modes do not exist in the liter- 
ature, nor is it entirely clear what the relevant damping 
mechanism is. We therefore advocate considerable more 
research in this area, associated with developing an un- 
derstanding of the relevant damping mechanism of these 
modes as well as targeting gravitational wave searches. 
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